#!/usr/bin/env Rscript

library("plyr")
library("foreach")
library("emIRT")
library("pscl")
library("ggplot2")
library("scales")


## #########
## Set Paths
## #########

p2root <- "../"
p2estimates <- paste(p2root, "data/RDATA/estimates/", sep = "")
p2figs <- paste(p2root, "doc/figs/", sep = "")

vAllFiles <- dir(p2estimates, recursive = TRUE)
vAllFiles <- setdiff(vAllFiles, "ideal/hou-112.rda")

## ########
## Function
## ########

getX1hat <- function(idx) {
    thisout <- lOut[[idx]]
    thistype <- thisout$output$type
    thisxhat <- switch(thistype,
                       'fastest' = {
                           thisout$output$output$means$x[, 1]
                       },
                       'fastestBS' = {
                           thisout$output$output$means$x[, 1]
                       },
                       'fastestMP' = {
                           thisout$output$output$means$x[, 1]
                       },
                       'fastestLOW' = {
                           thisout$output$output$means$x[, 1]
                       },
                       'wnominate' = {
                           thisout$output$output$legislators$coord1D
                       },
                       'ideal' = {
                           thisout$output$output$xbar[, 1]
                       },
                       'oc' = {
                           thisout$output$output$legislators$coord1D
                       },
                   {
                       stop(paste("unknown estimator:", thistype))
                   }
                       )
    return(list(thisxhat,
                thistype
                )
           )
}


## ################
## Read in All Data
## ################

lOut <- foreach(afile = vAllFiles,
                .combine = c
                ) %do% {
                    print(afile)
                    load(paste(p2estimates,
                               afile,
                               sep = ""
                               ),
                         verbose = TRUE
                         )
                    nRes <- length(lOut)
                    for (idx in 1:nRes) {
                        lOut[[idx]]$rc <- NULL
                        lOut[[idx]]$output$rc <- NULL
                    }
                    return(lOut)
                }




## ########################
## Create DF with Meta Data
## ########################

dfMeta <- ldply(lOut,
                function(x) {
                    return(data.frame(cham = x$cham,
                                      sess = x$sess,
                                      type = x$output$type,
                                      mc = x$mc
                                      )
                           )
                }
                )

head(dfMeta)

dfMeta$idx <- 1:nrow(dfMeta)

dfMeta2 <- subset(dfMeta, mc == 1 & sess == 112)

dfMeta2$idx

dfPoints <- ldply(dfMeta2$idx,
                  function(idx) {
                      data.frame(idx = idx,
                                 type = lOut[[idx]]$output$type,
                                 sess = lOut[[idx]]$sess,
                                 cham = lOut[[idx]]$cham,
                                 x = getX1hat(idx)[[1]],
                                 lidx = 1:length(getX1hat(idx)[[1]])
                                 )
                  }
                  )

dfF <- subset(dfPoints, type == "fastest")

dfO <- subset(dfPoints, type %in% c("ideal", "oc", "wnominate"))

dfAll <- merge(x = dfO,
               y = dfF,
               by = c("sess", "cham", "lidx"),
               suffixes = c(".other", ".vi")
               )

head(dfAll)


## ########
## Leg Data
## ########

dfLeg <- subset(dfMeta, type == "wnominate" & mc == 1)

dfLeg2 <- ddply(dfLeg,
                c("idx"),
                function(X) {
                    return(data.frame(idx = X$idx,
                                      cham = X$cham,
                                      sess = X$sess,
                                      lidx = 1:nrow(lOut[[X$idx]]$output$output$legislators),
                                      names = gsub(" \\(.*", "", rownames(lOut[[X$idx]]$output$output$legislators)),
                                      party = gsub("\\s.*", "", gsub(".*\\(", "", rownames(lOut[[X$idx]]$output$output$legislators)))
                                      )
                           )
                }
                )

head(dfLeg2)


### ##############################
### Ad Hoc Correlation Calculation
### ##############################

dfHleg <- subset(dfLeg2, cham == "hou" & sess == 112)
isD <- dfHleg$party == "D"
isR <- dfHleg$party == "R"

dfHF <- subset(dfPoints,
               sess == 112 &
               type == "fastest" &
               cham == "hou"
               )

dfHW <- subset(dfPoints,
               sess == 112 &
               type == "wnominate" &
               cham == "hou"
               )

dfHI <- subset(dfPoints,
               sess == 112 &
               type == "ideal" &
               cham == "hou"
               )

dfHO <- subset(dfPoints,
               sess == 112 &
               type == "oc" &
               cham == "hou"
               )


abs(cor(dfHF$x[isD], dfHW$x[isD]))
abs(cor(dfHF$x[isD], dfHI$x[isD]))

abs(cor(dfHW$x[isD], dfHI$x[isD]))

abs(cor(dfHF$x[isD], dfHO$x[isD], method = "spearman"))
abs(cor(dfHW$x[isD], dfHO$x[isD], method = "spearman"))
abs(cor(dfHI$x[isD], dfHO$x[isD], method = "spearman"))



abs(cor(dfHF$x[isR], dfHW$x[isR]))
abs(cor(dfHF$x[isR], dfHI$x[isR]))

abs(cor(dfHW$x[isR], dfHI$x[isR]))

abs(cor(dfHF$x[isR], dfHO$x[isR], method = "spearman"))
abs(cor(dfHW$x[isR], dfHO$x[isR], method = "spearman"))
abs(cor(dfHI$x[isR], dfHO$x[isR], method = "spearman"))




dfSleg <- subset(dfLeg2, cham == "sen" & sess == 112)
isD2 <- dfSleg$party == "D"
isR2 <- dfSleg$party == "R"

dfSF <- subset(dfPoints,
               sess == 112 &
               type == "fastest" &
               cham == "sen"
               )

dfSW <- subset(dfPoints,
               sess == 112 &
               type == "wnominate" &
               cham == "sen"
               )

dfSI <- subset(dfPoints,
               sess == 112 &
               type == "ideal" &
               cham == "sen"
               )

dfSO <- subset(dfPoints,
               sess == 112 &
               type == "oc" &
               cham == "sen"
               )


abs(cor(dfSF$x[isD2], dfSW$x[isD2]))

abs(cor(dfSF$x[isD2], dfSI$x[isD2]))

abs(cor(dfSW$x[isD2], dfSI$x[isD2]))

abs(cor(dfSF$x[isD2], dfSO$x[isD2], method = "spearman"))
abs(cor(dfSW$x[isD2], dfSO$x[isD2], method = "spearman"))
abs(cor(dfSI$x[isD2], dfSO$x[isD2], method = "spearman"))



abs(cor(dfSF$x[isR2], dfSW$x[isR2]))
abs(cor(dfSF$x[isR2], dfSI$x[isR2]))

abs(cor(dfSW$x[isR2], dfSI$x[isR2]))

abs(cor(dfSF$x[isR2], dfSO$x[isR2], method = "spearman"))
abs(cor(dfSW$x[isR2], dfSO$x[isR2], method = "spearman"))
abs(cor(dfSI$x[isR2], dfSO$x[isR2], method = "spearman"))












## #####
## Plots
## #####

dfAll2 <- subset(dfAll, type.other != "oc")

head(dfAll2)

dfAll3 <- ddply(dfAll2,
                .variables = c("sess", "cham", "type.other"),
                function(df) {
                    return(cbind(df, rel = ifelse(cor(df$x.vi, df$x.other) > 0, 1, -1)))
                }
                )

dfAll3$x.other2 <- dfAll3$x.other * dfAll3$rel

dfAll4 <- merge(dfAll3, dfLeg2, by = c("cham", "sess", "lidx"))

dfAll4$type.other2 <- ifelse(dfAll4$type.other == "ideal",
                             "Bayesian MCMC",
                             "W-NOMINATE"
                             )

dfAll4$cham2 <- ifelse(dfAll4$cham == "hou", "House", "Senate")


dfAll5 <- ddply(dfAll4,
      .variables= c("sess", "cham", "type.other"),
      function (df) {
          df2 <- within(df,
                    {
                        x.other2R <- (x.other2 - mean(x.other2)) / sd(x.other2)
                        x.viR <- (x.vi - mean(x.vi)) / sd(x.vi)
                    }
                        )
          return(df2)
      }
      )




## dfSumm1 <- ddply(dfAll4,
##                  c("cham2", "sess", "party", "type.other2"),
##                  function(X) data.frame(cor = cor(X$x.vi, X$x.other2))
##                  )


dfSumm1 <- ddply(dfAll5,
                 c("cham2", "sess", "party", "type.other2"),
                 function(X) data.frame(cor = cor(X$x.viR, X$x.other2R))
                 )


dfSumm2 <- ddply(dfAll5,
                 c("cham2", "sess", "type.other2"),
                 function(X) data.frame(cor = cor(X$x.viR, X$x.other2R),
                                        party = "all"
                                        )
                 )

pdf(file = paste(p2figs, "points_base.pdf", sep = ""),
    width = 10,
    height = 10
    )

par(mfcol = c(2, 2))

for (t in c("ideal", "wnominate")) {
    for (ch in c("hou", "sen")) {
        ch2 <- ifelse(ch == "hou", "House", "Senate")
        t2 <- ifelse(t == "ideal", "IDEAL", "W-NOMINATE")
        t3 <- ifelse(t == "ideal", "Bayesian MCMC", "W-NOMINATE")
        thisdf <- subset(dfAll5, cham == ch & type.other == t & party != "Indep")
        thisdf2 <- subset(dfSumm1, cham2 == ch2 & party != "Indep" & type.other2 == t3)
        print(thisdf2)
        plot(x = thisdf$x.other2R,
             y = thisdf$x.viR,
             xlim = c(-3.5, 3.5),
             ylim = c(-3.5, 3.5),
	     asp = 1,
             main = paste(ch2, ": ", t2, sep = ""),
             xlab = paste("Estimate (", ifelse(t == "ideal", "IDEAL", "W-NOMINATE") ,")", sep = ""),
             ylab = "Estimate (EM)",
             cex = .5,
             pch = ifelse(thisdf$party == "D", 1, 3.5),
             col = ifelse(thisdf$party == "D", alpha(muted("blue", l = 60), .85), alpha(muted("red", l = 60), .85)),
             cex.axis = 1.1,
	     cex.lab = 1.2,
             xaxt= "n",
             yaxt = "n"
             )
        axis(1, at = c(-2, 0, 2))
        axis(2, at = c(-2, 0, 2))
        text(x = 0,
             y = -2.5,
             labels = (substitute(paste(a, ": ", rho, " = ", b),
                                  list(a = "Democrats",
                                       b = round(thisdf2$cor[1], 3)
                                       )
                                  )
                       ),
             cex = 1.0
             )
        text(x = 0,
             y = 2.5,
             labels = (substitute(paste(a, ": ", rho, " = ", b),
                                  list(a = "Republicans",
                                       b = round(thisdf2$cor[2], 3)
                                       )
                                  )
                       ),
             cex = 1.0
             )
    }
}

dev.off()

################################################################################
